#pragma once

using namespace System;
using namespace NUnit::Framework;
using namespace LatoolNet;

namespace LatoolNetTest {

	[TestFixture]
	public ref class DoubleSymmetricMatrixTest {
	public:

		[Test]
		void TestPositiveDefiniteInvert() {
			
			int rownum = 4;
			int colnum = 4;

			Matrix ^a = gcnew Matrix(rownum, colnum, MatrixType::DoubleSymmetric);
			a[0, 0] = 4.16;
			a[0, 1] = -3.12;
			a[0, 2] = 0.56;
			a[0, 3] = -0.10;

			a[1, 1] = 5.03;
			a[1, 2] = -0.83;
			a[1, 3] = 1.18;

			a[2, 2] = 0.76;
			a[2, 3] = 0.34;

			a[3, 3] = 1.18;

			Matrix ^ ainv = a->Clone()->Inv();

			Matrix ^ uni = a * ainv;

			for (int i = 0; i < rownum; i++) {
				for (int j = 0; j < colnum; j++) {
					if (i == j) {
						Assert::AreEqual(1.0000, uni[i, j], 1e-5, "Test: Positive-definite Invert.");
					} else {
						Assert::AreEqual(0.0000, uni[i, j], 1e-5, "Test: Positive-definite Invert.");
					}
				}
			}

			Matrix ^ ainvinv = ainv->Clone()->Inv();

			for (int i = 0; i < rownum; i++) {
				for (int j = 0; j < colnum; j++) {
					Assert::AreEqual(a[i, j], ainvinv[i, j], 1e-5, "Test: Positive-definite Invert.");
				}
			}


		};


		[Test]
		void TestPositiveDefiniteFactorizeAndSolve() {
			
			int rownum = 4;
			int colnum = 4;

			Matrix ^a = gcnew Matrix(rownum, colnum, MatrixType::DoubleSymmetric);
			a[0, 0] = 4.16;
			a[0, 1] = -3.12;
			a[0, 2] = 0.56;
			a[0, 3] = -0.10;

			a[1, 1] = 5.03;
			a[1, 2] = -0.83;
			a[1, 3] = 1.18;

			a[2, 2] = 0.76;
			a[2, 3] = 0.34;

			a[3, 3] = 1.18;

			Matrix ^ b = gcnew Matrix(rownum, 1);

			b[0, 0] = 8.70;
			b[1, 0] = -13.35;
			b[2, 0] = 1.89;
			b[3, 0] = -4.14;

			//LUFactorization::Factorize(a);
			LinearEquation::Factorize(a);

			//LUFactorization::Solve(a, b);
			LinearEquation::Solve(a, b);

			//Console::WriteLine(b->ToString());
			Assert::AreEqual(1.0000, b[0, 0], 1e-10, "Test: Positive-definite Factorize And Solve.");
			Assert::AreEqual(-1.0000, b[1, 0], 1e-10, "Test: Positive-definite Factorize And Solve.");
			Assert::AreEqual(2.0000, b[2, 0], 1e-10, "Test: Positive-definite Factorize And Solve.");
			Assert::AreEqual(-3.0000, b[3, 0], 1e-10, "Test: Positive-definite Factorize And Solve.");

		};

		[Test]
		void TestIndefiniteFactorizeAndSolve() {
			
			int rownum = 4;
			int colnum = 4;

			Matrix ^a = gcnew Matrix(rownum, colnum, MatrixType::DoubleSymmetric);
			
			a[0, 0] = -1.81;
			a[0, 1] = 2.06;
			a[0, 2] = 0.63;
			a[0, 3] = -1.15;

			a[1, 1] = 1.15;
			a[1, 2] = 1.87;
			a[1, 3] = 4.20;

			a[2, 2] = -0.21;
			a[2, 3] = 3.87;

			a[3, 3] = 2.07;

			Matrix ^ b = gcnew Matrix(rownum, 1);

			b[0, 0] = 0.96;
			b[1, 0] = 6.07;
			b[2, 0] = 8.38;
			b[3, 0] = 9.50;

			//LUFactorization::Factorize(a);
			LinearEquation::Factorize(a);

			//LUFactorization::Solve(a, b);
			LinearEquation::Solve(a, b);

			Assert::AreEqual(-5.0000, b[0, 0], 1e-10, "Test: Indefinite Factorize And Solve.");
			Assert::AreEqual(-2.0000, b[1, 0], 1e-10, "Test: Indefinite Factorize And Solve.");
			Assert::AreEqual(1.0000, b[2, 0], 1e-10, "Test: Indefinite Factorize And Solve.");
			Assert::AreEqual(4.0000, b[3, 0], 1e-10, "Test: Indefinite Factorize And Solve.");

		};

		[Test]
		void TestPositiveDefiniteSolve() {
			
			int rownum = 4;
			int colnum = 4;

			Matrix ^a = gcnew Matrix(rownum, colnum, MatrixType::DoubleSymmetric);
			a[0, 0] = 4.16;
			a[0, 1] = -3.12;
			a[0, 2] = 0.56;
			a[0, 3] = -0.10;

			a[1, 1] = 5.03;
			a[1, 2] = -0.83;
			a[1, 3] = 1.18;

			a[2, 2] = 0.76;
			a[2, 3] = 0.34;

			a[3, 3] = 1.18;

			Matrix ^ b = gcnew Matrix(rownum, 1);

			b[0, 0] = 8.70;
			b[1, 0] = -13.35;
			b[2, 0] = 1.89;
			b[3, 0] = -4.14;

			//LUFactorization::Solve(a, b);
			LinearEquation::Solve(a, b);

			Assert::AreEqual(1.0000, b[0, 0], 1e-10, "Test: Positive-definite Solve.");
			Assert::AreEqual(-1.0000, b[1, 0], 1e-10, "Test: Positive-definite Solve.");
			Assert::AreEqual(2.0000, b[2, 0], 1e-10, "Test: Positive-definite Solve.");
			Assert::AreEqual(-3.0000, b[3, 0], 1e-10, "Test: Positive-definite Solve.");

		};

				[Test]
		void TestIndefiniteInvert() {
			
			int rownum = 4;
			int colnum = 4;

			Matrix ^a = gcnew Matrix(rownum, colnum, MatrixType::DoubleSymmetric);
			
			a[0, 0] = -1.81;
			a[0, 1] = 2.06;
			a[0, 2] = 0.63;
			a[0, 3] = -1.15;

			a[1, 1] = 1.15;
			a[1, 2] = 1.87;
			a[1, 3] = 4.20;

			a[2, 2] = -0.21;
			a[2, 3] = 3.87;

			a[3, 3] = 2.07;

			Matrix ^ ainv = a->Clone()->Inv();

			Matrix ^ ainvinv = ainv->Clone()->Inv();

			for (int i = 0; i < rownum; i++) {
				for (int j = 0; j < colnum; j++) {
					Assert::AreEqual(a[i, j], ainvinv[i, j], 1e-5, "Test: Indefinite Invert.");
				}
			}

			//Matrix ^ uni = a * ainv;

			//for (int i = 0; i < rownum; i++) {
			//	for (int j = 0; j < colnum; j++) {
			//		if (i == j) {
			//			Assert::AreEqual(1.0000, uni[i, j], 1e-5, "Test: Indefinite Invert.");
			//		} else {
			//			Assert::AreEqual(0.0000, uni[i, j], 1e-5, "Test: Indefinite Invert.");
			//		}
			//	}
			//}

		};


		[Test]
		void TestIndefiniteSolve() {
			
			int rownum = 4;
			int colnum = 4;

			Matrix ^a = gcnew Matrix(rownum, colnum, MatrixType::DoubleSymmetric);
			
			a[0, 0] = -1.81;
			a[0, 1] = 2.06;
			a[0, 2] = 0.63;
			a[0, 3] = -1.15;

			a[1, 1] = 1.15;
			a[1, 2] = 1.87;
			a[1, 3] = 4.20;

			a[2, 2] = -0.21;
			a[2, 3] = 3.87;

			a[3, 3] = 2.07;

			Matrix ^ b = gcnew Matrix(rownum, 1);

			b[0, 0] = 0.96;
			b[1, 0] = 6.07;
			b[2, 0] = 8.38;
			b[3, 0] = 9.50;

			//LUFactorization::Solve(a, b);
			LinearEquation::Solve(a, b);

			Assert::AreEqual(-5.0000, b[0, 0], 1e-10, "Test: Indefinite Solve.");
			Assert::AreEqual(-2.0000, b[1, 0], 1e-10, "Test: Indefinite Solve.");
			Assert::AreEqual(1.0000, b[2, 0], 1e-10, "Test: Indefinite Solve.");
			Assert::AreEqual(4.0000, b[3, 0], 1e-10, "Test: Indefinite Solve.");

		};

	};
}
